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The  Shock  Wave  Produced  by  Collapse 
of  a  Spherical  Cavity 
by 

R,  S.  Brand 

ABSTRACT 

The  results  of  the  computations  for  a  single  case  of  spherical 
cavity  collapse  in  water  are  presented,  including  the  post-collapse  shock 
wave  and  accompanying  velocity  and  pressure  fields.  For  the  model  studied 
it  appears  that  at  a  finite  time  after  collapse  negative  pressures  occur, 
indicating  that  the  cavity  reopens. 
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The  Shock  Wave  Produced  by  Collapse 
of  a  Spherical  Cavity 


PURPOSE  AND  SCOPE 


This  investigation  is  concerned  with  the  pressure  and  velocity  fields 
which  surround  an  empty  cavity  collapsing  with  spherical  symmetry  in  water. 

It  is  assumed  that  the  cavity  suddenly  appears  in  an  infinite  expanse 
of  water  which  is  at  rest  and  at  uniform  pressure,  p^, ,  Since  all  lengths 
are  made  dimensionless  with  respect  to  the  initial  cavity  radius,  the  only 
parameter  which  affects  the  collapse  is  ,  In  this  report  only  a  single 
case  is  presented,  that  for  which  p  .r  /B  =  0.02.  B  is  a  constant  appearing 
in  the  equation  of  state;  for  water  it  is  approximately  3  kilobars. 

The  water  is  assumed  to  be  inviscid  and  compressible. 

The  discontinuity  in  pressure  at  the  cavity  boundary  (zero  pressure 
:m  the  cavity,  p  just  outside)  cannot  persist.  An  expansion  wave  moves 
outward,  and  behind  this  wave  the  fluid  and  the  cavity  wall  move  inward. 

The  method  of  characteristics  is  used  to  compute  this  motion. 

Since  the  numerical  solution  cannot  be  carried  to  zero  radius,  it 
is  assumed  that  the  cavity  collapses  on  a  small  solid  sphere.  In  the  case 
''•’ported  here  the  ratio  of  the  radius  of  this  sphere  to  the  initial  cavity 
. s  , CO/*? 3 ,  When  the  inrushing  water  strikes  this  small  rigid  sphere,  its 
motion  is  arrested  by  a  spherical  shock  wave  moving  outward.  It  is  the 
major  object  of  this  investigation  to  describe  this  shock  wave.  It  is 
speculated  that  the  properties  of  the  shock  do  not  depend  strongly  on  the 
size  of  the  solid  sphere  once  the  shock  has  moved  a  small  distance  away. 


In  fact,  it  is  hoped  that  the  present  results  are  at  least  qualitatively 
descriptive  of  the  chock  that  would  be  produced  by  an  absolutely  empty 
cavity. 

Calculation  of  the  shock  is  carried  to  a  time  when  the  radius  of  the 
shock  is  almost  twice  the  initial  cavity  radius.  By  this  time  the  velocity 
of  the  shock  is  nearly  sonic  and  the  pressure  ratio  across  it  is  nearly 
one.  That  is,  it  has,  within  the  limits  of  accuracy  of  the  numerical  work, 
degenerated  into  a  continuous  wave. 

The  flow  behind  the  shock  is  also  computed  by  the  method  of  character¬ 
istics.  It  is  found  that  at  the  solid  sphere  the  pressure  eventually  falls 
to  zero  indicating  that  the  cavity  reopens.  Indeed,  a  sequence  of  several 
openings  and  closings  is  indicated  by  the  calculations.  However,  since 
the  principal  purpose  was  to  compute  the  shock,  the  behavior  at  the  solid 
boundary  has  not  been  adequately  investigated. 
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constant  appearing  in  the  equation  of  state 
speed  of  sound 

speed  of  sound  corresponding  to  zero  pressure 
dimensionless  speed  of  sound 
speed  of  sound  in  undisturbed  water 
dimensionless  speed  of  sound  just  outside  shock 
dimensionless  speed  of  sound  just  inside  shock 
pressure 
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dimensionless  pressure  in  undisturbed  water 
radial  distance  from  center  of  cavity- 
initial  cavity  radius 
dimensionless  radial  co-ordinate 
dimensionless  radius  of  solid  sphere 
time 

dimensionless  time 
velocity  of  fluid 
dimensionless  velocity  of  fluid 
velocity  of  shock  wave 
dimensionless  velocity  of  shock  wave 
exponent  in  equation  of  state 
density 

density  corresponding  to  zero  pressure 


THEORY 


It  is  assumed  that  at  time  t=€  an  empty  spherical  cavity  of  radius 
r^  suddenly  appears  in  an  infinite  expanse  of  water  which  is  at  rest  and 
at  some  uniform  pressure  p  ^  „ 

The  water  is  taken  to  be  a  compressible  fluid  having  an  equation  of 

state 


Y 


where  is  the  density  corresponding  to  pressure  p,  j->Q  is  the 
non-zero  density  corresponding  to  zero  pressure,  and  B  and  V  are  constants 
characteristic  of  the  liquid.  For  water  B  is  approximately  3  kilobars  and  V 
is  approximately  7„ 

The  speed  of  sound,  c,  can  be  computed  from  (l)  as 

_  sjifJL')  -  z  f_LY~' 

-  "fc  l  fj  '  l  ft  1  (2) 

where  co  is  the  sound  speed  that  would  prevail  in  water  at  zero  pressure. 
Hence  (1)  can  be  written 

P/B  -  (c/cQ)  2Y/(r  -1}  -1  (3) 


in  which  P  -  p/b  and  G  “  c/cq. 

The  dimensionless  sound  speed,  C,  is  the  most  convenient  thermodynamic 
state  variable,  and  all  results  are  given  in  terms  of  C.  The  pressure  is 
easily  computed  from  C  by  means  of  equation  (3). 
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It  is  apparent  that  relatively  small  changes  in  C  denote  very  high 
pressures.  For  example  C  =  1,03  corresponds  approximately  to  P  =  1.07 
or  p  =  210  bars.  (Recall  that  C  =  1,00  corresponds  to  p  »  0.)  In  some 
parts  of  the  computations  values  of  C  of  50  and  greater  were  obtained.  It 
is  doubtful  that  the  assumed  equation  of  state  is  valid  for  such  conditions. 

The  problem  is  to  determine  u,  the  fluid  velocity  and  c,  the  sound 
speed  as  functions  of  r,  the  distance  from  the  center  of  the  cavity,  and 
t,  the  time  measured  from  the  start  of  the  collapse.  These  variables  are 
made  dimensionless  by  choosing  cq  as  a  reference  velocity  and  rQ,  the  initial 
cavity  radius,  as  a  length  scale. 

R  -  r/rQ  T  =  cQt/ro  U  =  u/cq  C  =  c/cq  (4) 

In  terms  of  these  variables  the  differential  equations  expressing 
the  conservation  of  mass  and  of  momentum  can  be  combined  to  form  a  pair  of 
simultaneous  first  order  partial  differential  equations  in  which,  in  each 
equation,  the  dependent  variables  U  and  C  are  differentiated  in  the  same 


direction  in  the  R-T  plane. 


Q  .  -L-  [K  .. 

'  dR  Y-l  [dT  [  dR 


R 


* 


(7) 
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while  for  an  ingoing  characteristic  for  which 
dR  =  (U  -  C)  dT 

d“  -  A  dC  -  it  dT 


U  -  C,  we  have  from  (6) 


(8) 


The  computations  reported  here  are  based  on  finite  difference 

approximations  to  equations  (7)  and  (8). 

It  is  of  course  not  possible  to  carry  the  computations  all  the  way  to 

R  ~  0  because  of  the  ^  factors  occurring  on  the  right  hand  sides  of  (7)  and 

(8).  Accordingly  the  expedient  was  adopted  of  assuming  that  at  the  center 

of  the  cavity  there  exists  a  small  solid  sphere  of  (dimensionless)  radius 

R  R  .  was  not  chosen  at  the  outset  but  was  taken  to  be  the  smallest 

min  min 

positive  value  of  R  which  could  be  reached  by  reducing  the  mesh  size  of  the 
characteristic  network  a  pre-selected  number  of  times.  In  the  case  reported 
here  R^^  =  0.00423. 

When  the  inrushing  cavity  boundary  strikes  this  solid  sphere,  its 
motion  is  arrested  instantly.  This  discontinuity  in  the  fluid  velocity  is 
the  origin  of  a  shock  wave  which  moves  outward  very  rapidly  at  first  but 
with  decreasing  velocity.  The  shock  wave  is  treated  here  as  a  true 
mathematical  discontinuity;  that  is,  at  the  shock  at  each  instant  there  are 
two  values  of  U  and  two  of  C.  These  are  related  to  each  other  and  to  W 


the  dimensionless  shock  velocity  by 

~  ^  _  /  CaY^ 

uA-w  \cBJ 


(U8-w)(uA-u8)  =  -i-c* 


(9) 

(10) 
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in  which  subscript  11  A"  refers  to  the  fluid  just  ahead  (i.e.  outside)  the 
shock  and  subscript  “B"  refers  to  the  fluid  inside  the  shock. 

Equations  (9)  and  (10)  are  derived  from  the  statements  of  conservation 
of  mass  and  momentum  together  with  the  assumption  that  the  equation  of  state 

(l)  holds  across  the  shock,  an  assumption  justified  by  the  fact  that  the  state 
of  water  is  quite  insensitive  to  entropy.  Hence  the  mechanical  conditions 
are  sufficient  to  describe  the  shock  process  and  no  consideration  need  be 
given  to  the  energy  equation.  (See  reference  1,  pp.  131-132) 

Solutions  to  equations  (9)  and  (10)  must  be  obtained  at  every  point 
on  the  shock  which  are  compatible  with  the  characteristic  equations  on 
both  sides.  A  more  detailed  description  of  how  this  was  done  is  given  in 
the  Appendix. 

The  boundary  conditions  for  the  problem  are  as  follows: 

(1)  On  R  =  T  +  1.  which  is  the  leading  edge  of  the  outgoing  expansion  wave 

which  originates  at  R  =  1.0,  T  =  0,  we  have  P  =  or  0  =  C^,  . 

(2)  During  the  collapse  phase  on  the  cavity  boundary,  which  is  an 
undetermined  curve  in  the  R-T  plane,  C  =  1 

(3)  After  collapse,  when  R  ■=  R^^,  U  =  0 

It  was  found  that  in  applying  boundary  condition  (3),  values  of  C<  1.0 
were  encountered.  This  would  imply  that  the  cavity  was  reopened.  When 
this  occurred,  the  values  obtained  were  rejected  and  boundary  condition  (2) 
was  substituted  for  (3). 

Several  cycles  of  opening  and  closing  were  obtained  but  the  size  of 
the  finite  steps  taken  was  such  that  this  behavior  is  not  well  delineated. 


RESULTS 


The  results  of  the  computations  are  shown  in  the  curves.  Figs,  (l) 
through  (8). 

Fig.  (l)  shows  the  cavity  wall  path,  the  shock,  and  several 
representative  characteristics.  Figs.  (2)  and  (3)  are  enlargements  near 
the  collapse  point  of  the  same  information.  It  should  be  pointed  out  that 
even  in  Fig.  (3),  which  shows  considerable  detail  near  R  *=  Rm^n,  only  a 
fraction  of  the  characteristics  actually  computed  have  been  plotted. 

Fig.  (4)  is  a  log-log  plot  of  the  velocity  of  the  cavity  boundary 
vs.  the  cavity  radius.  From  the  slope  of  this  curve,  which  is  nearly  a 
straight  line  for  small  R,  it  can  be  determined  that  the  final  collapse 
appears  to  follow  a  law  of  the  type. 


with  d,  =  -0.90.  This  tends  to  confirm  Hunter's  analysis  (Reference  2) 

Figs.  (5)  and  (6)  are  descriptive  of  the  shock,  showing  its  velocity 
and  strength  (i.e.  the  sound  speed  on  both  sides  of  the  shock)  plotted 
against  the  radius  of  the  shock. 

Figs.  (7)  and  (8)  show  the  variation  of  sound  speed  with  radius  at 
different  times,  and  is  thus  descriptive  of  the  pressure  field. 

All  of  the  results  given  in  the  curves  are  for  =  p,„ /B  =  0.02 

and  R  .  -  0.00423 
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APPENDIX 


Details  of  the  Computations 


A,  Characteristics 

Computation  proceeded  stepwise  along  the  ingoing  characteristics. 
Assume  that,  in  Fig.  (9),  values  of  H,  T,  U,  and  C,  have  been  found  at 
all  the  mesh  points  indicated  by  the  solid  lines,  including  points  1  and 
2,  and  we  now  wish  to  compute  R^,  T^,  U^,  and  C^. 

The  finite  difference  approximations  to  (7)  and  (8)  are 

r3  -  \  =  J(ux  +  c1  +  u3  +  c3)  (t3  -  Tj)  (11) 

U3  "  U1  +  yII  (C3  ~  Ci)  “  “  (U1C1  +  (T-  -  Tx)  (12) 

R1  R3 

R3  -  R2  =  £( U2  ~  c2  +  u3  -  c3)  (t3  -  T2)  (13) 

U3  "  U2  ~  yli  "  C2)  “  (U2°2  +  ^&)  (T3  -  T2)  (14) 

R~  R- 


The  simultaneous  solution  of  these  four  non-linear  equations  was 
effected  by  an  iterative  process  on  the  equations 


R3  -  R, 


K1  (T3  "  V 


U3  “  U1  4  f-1  (C3  -  Cl} 


K2  (T  -  Tx) 


r3-r2 


k3  (t3  -  t2) 


u3  -  u2  -  7H  (C3 


w,  -  CJ  =  K4  (T3  -  T J 


(15) 

(16) 

(17) 

(13) 


Initially  the  coefficients  were  approximated  by 


K1  -  Ul  +  Cl 


U  G 

K  =  2  — 

2  R, 


u2-c2 


U2C2 

VHt 


all  of  which  are  Known.  Equations  (15)  through  (18)  are  then  linear  and 
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are  easily  solved  for  R^,  T^,  U^,  C^.  New  values  of  the  coefficients 
are  then  computed  from 

u.c  uc 

K!  "  ^U1  +  C1  +  U3  +  C3)  K2  ”  rT  +  rT 


k3  -  i(u2  -  C2  +  u3  -  c3) 


K  +  ^ 

2  R1 

K  + 

4  r2  r3 


Use  of  these  improved  values  of  the  in  equations  (15)  through  (18) 
lead  to  improved  values  of  R^,  T C y  The  process  is  repeated  until 
the  change  in  R3  from  one  iteration  to  the  next  is  less  than  a  pre-selected 
tolerance. 


B.  The  Shock  Wave 


Let  it  be  assumed  that  calculations  have  been  completed  on  the 
characteristic  labelled  "100"  in  Fig.  (10) ,  including  all  values  at  the 
shock.  Calculations  on  characteristic  ”101“  have  been  carried  up  to  and 
Including  point  U,  and  wo  now  wish  to  compute  the  location  of  and  flow 
parameters  at  point  5.  That  is,  the  values  of  R,  T,  U,  and  C  are  known 
at  points  1,  3,  and  4,  and  at  2  we  have  R,  T,  U^,  C^,  Ug,  Cg,  and  W.  We 
now  wish  to  calculate  the  seven  unknowns  associated  with  the  shock  at 


point  5  o 

Point  5  is  at  the  intersection  of  the  shock  wave  with  the  ingoing 
characteristic  “101".  Because  a  shock  wave  always  has  a.  supersonic  velocity 
relative  to  the  fluid  ahead  of  it,  there  must  be  an  outgoing  characteristic 
arising  from  some  point  on  characteristic  “lOO"  between  points  1  and  2, 
which  also  passes  through  point  5.  This  auxiliary  point  is  labelled  7  in 
the  diagram  and  the  auxiliary  characteristic  is  shown  as  a  dashed  line. 

Similarly,  behind  the  shock,  there  must  be  a  point  8,  between  2  and 
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3,  where  an  outgoing  characteristic  will  also  intersect  the  shock  at  point 
5S  because  a  shock  always  has  a  subsonic  speed  relative  to  the  fluid  behind  it. 

It  is  assumed  that  the  parameters  at  point  7  are  related  by  a  linear 

interpolation  between  points  1  and  2,  and  that  8  is  similarly  related  to  2  and  3. 

We  now  summarize  the  unknowns  and  the  available  equations. 

Unknowns i 

At  point  5  -  r5»  t5»  ua_5>  cA_5>  UB-5>  CB-5j  W5 

At  point  7  -  R^,  T^,  U^, 

At  point  8  -  Rg(  Tg,  Ug,  Cg 

Equations 

Shock  relations  (9)  and  (10)  at  point  5  (2  equations) 

Characteristic  equations  relating  R^,  T^,  UA  ^  to 

points  4  and  7.  (4  equations) 

Characteristic  equations  for  outgoing  characteristic  only,  relating 

equations) 

Interpolation  for  point  7  between  1  and  2  (3  equations) 

Interpolation  for  point  8  between  2  and  3  (3  equations) 

Velocity  of  the  shock 

(R5  -  R2)  /  (T5  -  T2)  -  J(W5  +  W2)  (1  equation) 

Thus  we  have  15  relationships  available  for  determing  the  15 
unknowns „ 

The  sequence  adopted  for  solving  these  equations  was  as  follows: 

1 .  Assume  a  value  of  (The  initial  assumption  was  W^  =  W,,) 

2.  Assume  a  value  of  R,,,  (Initial  choice:  R^  ■=  £(R^  +  R2) 

3.  Find  T^j,  ,  and  C^  by  interpolation  between  1  and  2. 
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4.  Using  points  4  and  7  as  input,  find  R^,  T^,  UA  <.,  by  means  of  the 

characteristics  subroutine. 

5.  Test  whether  the  point  R^,  T,.  determined  in  step  4  lies  on  the  shock; 

(»5  -  r2)  /  (t5  -  t2)  -  i(w2  +  W$)  ? 

If  this  relation  is  not  satisfied  closely  enough,  adjust  the  value  of 
R^  and  return  to  step  3» 

If  the  relation  is  satisfied,  proceed  to  step  6. 

6.  Since  we  now  have  values  of  U,  CA  r  and  Wr,  the  shock  relations  have 

A-5S  A-5  5* 

only  two  unknowns,  Ug  and  Cg  ^ „  Solve  these  simultaneously.  (A  Newton- 
Raphson  iteration  was  used,  rather  than  simple  elimination  of  one  of  the 
unknowns  because  of  the  complexity  of  the  equations.) 

7.  Assume  a  value  of  Rg  (Initial  choice:  Rg  =  £(R2  +  ^3) 

8.  Interpolate  for  Tg,  Ug,  and  Cg,  between  points  2  and  3 

9.  Test  whether  the  point  R^,  T^,  lies  on  the  outgoing  characteristic  through 
point  8; 

(R5  -  Rs)  /  <t5  -  T8)  -  k%  *  c8  *  uB_5  *  cB_5)  ? 

If  this  test  is  not  met  within  the  specified  tolerance,  adjust  R0  and 

O 

return  to  step  8. 

If  the  test  is  met,  proceed  to  step  10, 

10.  Compare  the  value  of  determined  by  the  shock  in  step  6  with  the 

value  of  Ug  _  that  the  characteristic  relation  for  am  outgoing  character¬ 
istic  would  predict.  Let  the  difference  between  these  values  of  UD  c  be 
an  error  E 


E  =  U 


“a  *  A  lcB-5  -  V  -  <*5  -  V 


B-5  (shock) 


11.  Assume  a  new  value  of  somewhat  less  than  and  repeat  the  entire 
process  of  steps  2  through  10,  finding  a  new  value  of  E  corresponding 
to  the  new  value  of 

12.  Assuming  a  linear  relationship  between  E  and  W,.,  make  a  third  choice 

of  W.  which  will  reduce  E  to  zero. 

5 

13.  Continue  until  the  (n  +  l)^*1  choice  of  differs  from  the  nth  choice  by 
less  than  a  pre-assigned  tolerance. 

C.  Points  on  the  Cavity  Boundary 

On  the  boundary  of  the  collapsing  cavity  the  pressure  is  zero  and 
hence  C  =  1.  The  velocity  of  the  liquid  i3  the  velocity  of  the  boundary  itself. 

In  Fig.  11,  let  us  assume  that  the  calculations  for  ingoing  character¬ 
istic  "50'*  are  complete  (including  the  boundary  point,  l),  and  that  on 
characteristic  '*51"  we  have  computed  up  to  and  including  point  2  by  the  method 
of  characteristics  as  described  in  Section  A  above.  We  now  wish  to  establish 
boundary  point  3  on  characteristic  ”51". 

Unknowns  are  R^,  T^,  U^.  is  known  to  be  unity  because  of  the 
condition  of  zero  pressure  in  the  cavity. 

The  fact  that  both  points  1  and  3  are  on  the  cavity  boundary  where 
the  velocity  of  the  boundary  is  the  velocity  of  the  fluid  is  expressed  by 
(R3  -  R1)  /  (T3  -  Tx)  =  +  U3)  (19) 

Since  both  3  and  2  are  on  the  same  ingoing  characteristic,  relations 
(13)  and  (ll)  apply 

R3  -  R2  -  £(U2  -  c2  +  u3  -  c3)  (t3  -  t2)  (20) 

2  u  c  U  C 

u3  -  u2  -  y?i  (c3  -  C2}  -  (-tr +  v  (t3  -  v  (21) 

Equations  (19) ,  (20)  and  (21)  are  sufficient  to  determine  R3>  7^t 
and  U  .  Actual  solution  of  the  equations  was  carried  out  by  an  iterative 


technique  very  similar  to  that  employed  in  the  method  of  characteristics  for 
points  not  on  the  boundary. 

D.  Points  on  the  Solid  Boundary 
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The  calculations  are  based  on  the  assumption  that  the  cavity  collapses 


onto  a  rigid  sphere  of  radius  Rm^n.  Hence  the  ingoing  characteristics  which 
cross  the  shock  are  terminated  at  that  radius,  and  in  Fig.  12  it  is  necessary 
to  compute  point  2  from  point  1,  using  R2  "  Rmin>  U2  =  0}  with  Cg  and 
determined  by  the  ingoing  characteristic  equations: 

R2  -  R1  ’  i<U2  *  C2  +  U1  "  V  <T2  "  V 


u2  -  U 


-  U-C,  lie. 

h  iC2  -  Cl>  -  (-kf  +  ^  ^  ^ 


Y-l  ""2  “r  '  R„  ■  R,  ‘  2  *1' 

Simultaneous  solution  of  these  equations  for  and  is  routine. 
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